Wind energy

Intro to wind energy

Exercise 1

Map wind turbines in Rogaland from OSM and calculate the total energy capacity installed.

Objetives

• OSM (and how get data from it)
• Map points data (interactive map)
  
• Manipulate spatial data (e.g., intersection points - polygons)
• Types of data in a data frame
• Data wrangling

Solution

Load libraries

library(osmdata)
library(sf)
library(tmap)
library(tidyverse)

Get data from OpenStreetMap with osmdata, and transform to a sf object.

query <- opq(bbox = "Rogaland", timeout = 50) |> 
   add_osm_feature(key = "generator:source", value = "wind")
data <-  osmdata_sf(query)

wind_turbines <- data$osm_points 

Names of the data frame,

wind_turbines |>
  names()

 [1] "osm_id"                           "name"                            
 [3] "aeroway:light"                    "aeroway:light:character"         
 [5] "aeroway:light:colour"             "aeroway:light:icao_type"         
 [7] "aeroway:light:intensity"          "brand"                           
 [9] "date_constructed"                 "description"                     
[11] "ele"                              "generator:method"                
[13] "generator:output:electricity"     "generator:source"                
[15] "generator:type"                   "height"                          
[17] "man_made"                         "model"                           
[19] "operator"                         "power"                           
[21] "ref"                              "ref:hinder"                      
[23] "ref:kystverket"                   "seamark:landmark:category"       
[25] "seamark:landmark:conspicuity"     "seamark:landmark:vertical_length"
[27] "seamark:light:1:character"        "seamark:light:1:colour"          
[29] "seamark:light:1:group"            "seamark:light:1:height"          
[31] "seamark:light:1:multiple"         "seamark:light:1:period"          
[33] "seamark:light:1:range"            "seamark:light:2:category"        
[35] "seamark:light:2:character"        "seamark:light:2:colour"          
[37] "seamark:name"                     "seamark:type"                    
[39] "serial_number"                    "source"                          
[41] "geometry"                        

Clean names

wind_turbines <- wind_turbines |>
  janitor::clean_names()

names(wind_turbines)

 [1] "osm_id"                           "name"                            
 [3] "aeroway_light"                    "aeroway_light_character"         
 [5] "aeroway_light_colour"             "aeroway_light_icao_type"         
 [7] "aeroway_light_intensity"          "brand"                           
 [9] "date_constructed"                 "description"                     
[11] "ele"                              "generator_method"                
[13] "generator_output_electricity"     "generator_source"                
[15] "generator_type"                   "height"                          
[17] "man_made"                         "model"                           
[19] "operator"                         "power"                           
[21] "ref"                              "ref_hinder"                      
[23] "ref_kystverket"                   "seamark_landmark_category"       
[25] "seamark_landmark_conspicuity"     "seamark_landmark_vertical_length"
[27] "seamark_light_1_character"        "seamark_light_1_colour"          
[29] "seamark_light_1_group"            "seamark_light_1_height"          
[31] "seamark_light_1_multiple"         "seamark_light_1_period"          
[33] "seamark_light_1_range"            "seamark_light_2_category"        
[35] "seamark_light_2_character"        "seamark_light_2_colour"          
[37] "seamark_name"                     "seamark_type"                    
[39] "serial_number"                    "source"                          
[41] "geometry"                        

Interactive map with tmap.

tmap_mode("view")

tm_shape(wind_turbines) +
  tm_dots(col = "#0072B2")

Figure 1: Wind turbines in Rogaland (Data from OSM)

There are some wind turbines that ar not in Rogaland, so we are going to delect them from the dataset. For that we need to intersect our points (wind turbines) with the polygon (Rogaland).

# Get Norwegian counties (polygons)) from GISCO
counties <- giscoR::gisco_get_nuts(country = "NO",
                                   year = "2021",
                                   nuts_level = 3,
                                   epsg = "4326",
                                   resolution = "01") 
# Get only rogaland county
rogaland <- counties |> 
  filter(NUTS_NAME == "Rogaland")

# Intersect wind turbines (points) with rogaland (polygons)
wind_turbines <- wind_turbines |> 
  st_intersection(rogaland)

# Plot
tm_shape(wind_turbines) +
  tm_dots(col = "#0072B2")

Figure 2: Wind turbines in Rogaland (Data from OSM)

Now, we are going to calculate the maximum capacity (MW) installed in the region. If we inspect the data, there is a column describing called generator:output:electricity, which give the information we are looking for in MW. However, if we inspect the columns we see that the values are characters (<chr>) not numbers, so we can not carried out numerical operations on them.

glimpse(wind_turbines)

Rows: 256
Columns: 50
$ osm_id                           <chr> "1277731149", "1510248163", "15102482…
$ name                             <chr> "Høg-Jæren energipark 31", "Høg-Jæren…
$ aeroway_light                    <chr> "obstacle", NA, "obstacle", "obstacle…
$ aeroway_light_character          <chr> "flashing", NA, "flashing", "flashing…
$ aeroway_light_colour             <chr> "red", NA, "red", "red", "red", "red"…
$ aeroway_light_icao_type          <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ aeroway_light_intensity          <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ brand                            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ date_constructed                 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ description                      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ ele                              <chr> "264", "268", "250", "234", "205", "2…
$ generator_method                 <chr> "wind_turbine", "wind_turbine", "wind…
$ generator_output_electricity     <chr> "2.3 MW", "2.3 MW", "2.3 MW", "2.3 MW…
$ generator_source                 <chr> "wind", "wind", "wind", "wind", "wind…
$ generator_type                   <chr> "horizontal_axis", "horizontal_axis",…
$ height                           <chr> "126", "126", "126", "126", "89", "89…
$ man_made                         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ model                            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ operator                         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ power                            <chr> "generator", "generator", "generator"…
$ ref                              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ ref_hinder                       <chr> "70647", "70646", "70640", "71472", "…
$ ref_kystverket                   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_landmark_category        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_landmark_conspicuity     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_landmark_vertical_length <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_character        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_colour           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_group            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_height           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_multiple         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_period           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_1_range            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_2_category         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_2_character        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_light_2_colour           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_name                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ seamark_type                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ serial_number                    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ source                           <chr> "Kartverket Luftfartshindre", "Kartve…
$ NUTS_ID                          <chr> "NO0A1", "NO0A1", "NO0A1", "NO0A1", "…
$ LEVL_CODE                        <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
$ URBN_TYPE                        <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CNTR_CODE                        <chr> "NO", "NO", "NO", "NO", "NO", "NO", "…
$ NAME_LATN                        <chr> "Rogaland", "Rogaland", "Rogaland", "…
$ NUTS_NAME                        <chr> "Rogaland", "Rogaland", "Rogaland", "…
$ MOUNT_TYPE                       <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ COAST_TYPE                       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ geo                              <chr> "NO0A1", "NO0A1", "NO0A1", "NO0A1", "…
$ geometry                         <POINT [°]> POINT (5.792002 58.65093), POIN…

Therefore, we need to transform the data to numbers. We can use the function parse_number() from the package readr (inside tidyverse). For not deleting the column, we can generate a new one (e.g., gen_electricity_mw)

wind_turbines <- wind_turbines |> 
  mutate(gen_electricity_mw = parse_number(generator_output_electricity) )
     
# Show first 10 values of the column
wind_turbines$gen_electricity_mw |> 
  head(10)

 [1] 2.3 2.3 2.3 2.3  NA  NA 2.3 2.3 2.3 2.3

Remove NA in energy generator:

wind_turbines <- wind_turbines |> 
  drop_na(gen_electricity_mw)

# Summary of energy generator (without NA) 
wind_turbines$gen_electricity_mw |> 
  head(10)

 [1] 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3 2.3

calculate the total generation energy capacity.

total_capacity_mw <- wind_turbines$gen_electricity_mw |> sum()

print(total_capacity_mw)

[1] 905.3

So the total wind capacity installed in Rogaland is 905.3 MW. Note, that this is only the maximum energy installed and not how much it is generate in the region. It does not take into account the efficiency of the turbines nor wind availability!!.

Exercise 2

Map wind farms and wind turbines from NVE (www.nve.no). The data are free but need to be downloaded from https://nedlasting.nve.no/gis/ before reading into R (save it in a folder: e.g., ~/data/big_data/NVE/NVEData). I have downloaded them in .geojson format. Therefore, we need to read them with the geojsonsf package, which converts GeoJSON to sf objects.

Objetives

• Load data from local files
• Data wrangling (preprocessing)
• Formats of spatial data
• Spatial intersections (i.e., points - polygons)
• Types of vector data (i.e., lines, points, polygons)
• Plot more that one layer in a interactive map

Solution

# Libraries
library(geojsonsf)
library(sf)
library(tidyverse)
library(tmap)

Load data from a local file.

1. Wind turbines (point data)

wind_turbines_nve_path <- "data/big_data/NVE/NVEData/Vindkraft_Vindturbin.geojson"
wind_turbines_nve <- geojson_sf(wind_turbines_nve_path) |> 
  # dataUttaksdato to date format
  mutate(dataUttaksdato = ymd(dataUttaksdato))

wind_turbines_nve

Simple feature collection with 272 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4.905257 ymin: 58.31454 xmax: 6.51592 ymax: 59.41634
Geodetic CRS:  WGS 84
First 10 features:
   objektType   saksTittel saksKategori dataUttaksdato              eksportType
1  Vindturbin      Svåheia            2     2024-02-07 NVEs nedlastningsløsning
2  Vindturbin       Tysvær            2     2024-02-07 NVEs nedlastningsløsning
3  Vindturbin    Høg-Jæren            2     2024-02-07 NVEs nedlastningsløsning
4  Vindturbin       Tysvær            2     2024-02-07 NVEs nedlastningsløsning
5  Vindturbin    Bjerkreim            2     2024-02-07 NVEs nedlastningsløsning
6  Vindturbin Vardafjellet            2     2024-02-07 NVEs nedlastningsløsning
7  Vindturbin    Bjerkreim            2     2024-02-07 NVEs nedlastningsløsning
8  Vindturbin   Måkaknuten            2     2024-02-07 NVEs nedlastningsløsning
9  Vindturbin       Utsira            2     2024-02-07 NVEs nedlastningsløsning
10 Vindturbin     Tellenes            2     2024-02-07 NVEs nedlastningsløsning
                    geometry
1  POINT (6.091838 58.38779)
2  POINT (5.552303 59.31521)
3   POINT (5.76918 58.64514)
4   POINT (5.557336 59.2931)
5  POINT (5.934883 58.59081)
6  POINT (5.894319 58.83959)
7  POINT (5.945595 58.58838)
8  POINT (5.966087 58.68458)
9  POINT (4.906018 59.31497)
10 POINT (6.514851 58.34872)

2. Wind farms areas (Polygons)

wind_farms_nve_path <- "data/big_data/NVE/NVEData/Vindkraft_VindkraftanleggOmr.geojson"
wind_farms_nve <- geojson_sf(wind_farms_nve_path) |> 
  # Coherce to dates format
  mutate(across(.cols = ends_with("dato"), .fns = ymd)) |> 
  # Add ID column
   rowid_to_column("id_farm")

wind_farms_nve

Simple feature collection with 57 features and 15 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 4.890143 ymin: 58.2376 xmax: 6.520988 ymax: 59.61306
Geodetic CRS:  WGS 84
First 10 features:
   id_farm dataUttaksdato totaltAntTurbiner    objektType utAvDriftDato
1        1     2024-02-07              <NA> Vindkraftverk          <NA>
2        2     2024-02-07              <NA> Vindkraftverk          <NA>
3        3     2024-02-07              <NA> Vindkraftverk          <NA>
4        4     2024-02-07              <NA> Vindkraftverk          <NA>
5        5     2024-02-07                37 Vindkraftverk          <NA>
6        6     2024-02-07              <NA> Vindkraftverk          <NA>
7        7     2024-02-07              <NA> Vindkraftverk          <NA>
8        8     2024-02-07              <NA> Vindkraftverk          <NA>
9        9     2024-02-07              <NA> Vindkraftverk          <NA>
10      10     2024-02-07              <NA> Vindkraftverk          <NA>
   idriftDato status fylkesnavn
1        <NA>      V   Rogaland
2        <NA>      V   Rogaland
3        <NA>      V   Rogaland
4        <NA>     FJ   Rogaland
5  2019-12-19      D   Rogaland
6        <NA>     FJ   Rogaland
7        <NA>     FJ   Rogaland
8        <NA>      V   Rogaland
9        <NA>      V   Rogaland
10       <NA>     FJ   Rogaland
                                                                    saksLenke
1  https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=194
2   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=29
3   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=36
4   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=52
5   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=32
6   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=37
7   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=64
8   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=55
9   https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=53
10 https://www.nve.no/konsesjon/konsesjonssaker/konsesjonssak?type=A-6&id=120
   forventetProduksjon_Gwh kommunenavn                       geometry
1                     34.0   Stavanger POLYGON ((5.616381 59.07174...
2                    408.0   Bjerkreim POLYGON ((5.892451 58.63726...
3                    785.4   Bjerkreim POLYGON ((5.893369 58.70034...
4                    272.0        Lund POLYGON ((6.26167 58.56276,...
5                    558.4   Bjerkreim POLYGON ((5.995785 58.58672...
6                    306.0   Eigersund POLYGON ((6.206857 58.56671...
7                    306.0     Sokndal POLYGON ((6.153781 58.36661...
8                    544.0     Gjesdal POLYGON ((6.293888 58.82249...
9                    265.2   Bjerkreim POLYGON ((6.07609 58.68459,...
10                   204.0      Karmøy POLYGON ((5.29149 59.23401,...
   saksKategori               tiltakshaver effekt_MW              eksportType
1             2 MARIN ENERGI TESTSENTER AS      10.0 NVEs nedlastningsløsning
2             3         LYSE PRODUKSJON AS     120.0 NVEs nedlastningsløsning
3             3         LYSE PRODUKSJON AS     231.0 NVEs nedlastningsløsning
4             4          BJERKREIM VIND AS      80.0 NVEs nedlastningsløsning
5             2          BJERKREIM VIND AS     159.1 NVEs nedlastningsløsning
6             4          BJERKREIM VIND AS      90.0 NVEs nedlastningsløsning
7             4                  ZEPHYR AS      90.0 NVEs nedlastningsløsning
8             2     GILJA VINDKRAFTVERK AS     160.0 NVEs nedlastningsløsning
9             3                  ZEPHYR AS      78.0 NVEs nedlastningsløsning
10            4      EQUINOR WIND POWER AS      60.0 NVEs nedlastningsløsning

Map both datasets together:

tmap_mode("view")


tm_shape(wind_farms_nve) + 
  tm_fill("status", alpha = 0.5, title = "Wind farm status") +
  # Add wind turbines
  tm_shape(wind_turbines_nve) +
  tm_dots(col = "#0072B2")

Figure 3: Wind turbines in Rogaland (Data from NVE)

Note: Status

D - Drift (Operations) N - Nedlagt (Decommissioned) O - Ombygd (Rebuilt) P - Planlagt (Planned) P1 - Planlagt illustrert (Planed illustrated) P2 - Planlagt, prosjektert (Planed, projected) U - Under arbeid (in progress) V - Vedtatt (Adopted) FJ - Fjernet (Removed)

To calculate the actual capacity installed, we need to select only the wind farms that are in operation from wind_farms_nve, and sum the power capacity (effekt_MW). We can do that by montds to see the temporal evolution.

power_year <- wind_farms_nve |>
  # Get only farm in operation
  filter(status == "D") |> 
  # Summarize power by month
  group_by(year = lubridate::floor_date(idriftDato, "year")) %>%
  summarize(power_MW = sum(effekt_MW)) |> 
  ungroup() |> 
  # Cummulative sum
  mutate(cumsum_power_MW = cumsum(power_MW))

# Column plot
ggplot(data = power_year,
       aes(x = year, y = cumsum_power_MW)) +
  geom_col(fill = "darkblue") +
  labs(title = "Cumulative wind power installed in Rogaland",
       y = "Power [MW]",
       x = "") +
  theme_bw()

Wind capacity

We can count now the number of wind turbines per wind farm, to understand for example the volume of blades we may need to recycler.

 # Intersect points (wind turbines) wit polygons (wind farms)
number_turbines_farm <- wind_turbines_nve |> 
  # Detect wind farm 
  st_intersection(wind_farms_nve) |> 
  # Number of turbines 
  group_by(id_farm, status,  idriftDato, effekt_MW) |> 
  summarize(n = n()) |> 
  ungroup()

number_turbines_farm

Simple feature collection with 28 features and 5 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 4.905257 ymin: 58.31454 xmax: 6.51592 ymax: 59.41634
Geodetic CRS:  WGS 84
# A tibble: 28 × 6
   id_farm status idriftDato effekt_MW     n                            geometry
     <int> <chr>  <date>         <dbl> <int>                      <GEOMETRY [°]>
 1       2 V      NA              120      1           POINT (5.885113 58.58739)
 2       3 V      NA              231     22 MULTIPOINT ((5.909491 58.66277), (…
 3       5 D      2019-12-19      159.    37 MULTIPOINT ((5.902099 58.58625), (…
 4      14 D      2019-07-25       90     18 MULTIPOINT ((5.871099 58.57491), (…
 5      16 D      2020-06-22       30      7 MULTIPOINT ((5.955698 58.66794), (…
 6      19 D      2020-08-17       30      7 MULTIPOINT ((5.88813 58.8301), (5.…
 7      29 D      2018-02-13       10      3 MULTIPOINT ((5.93595 58.72669), (5…
 8      30 D      2018-10-12       10      2 MULTIPOINT ((5.908438 58.74937), (…
 9      32 P2     NA               15      4 MULTIPOINT ((4.905952 59.31499), (…
10      33 D      2017-11-16      110     33 MULTIPOINT ((6.119956 58.4561), (6…
# ℹ 18 more rows

Now we are going to plot the evolution but in the number od witd turbines in operation:

turbines_year <- number_turbines_farm |>
  # Get only farm in operation
  filter(status == "D") |> 
  # Summarize power by month
  group_by(year = lubridate::floor_date(idriftDato, "year")) %>%
  summarize(num_turbines = sum(n)) |> 
  ungroup() |> 
  # Cummulative sum
  mutate(cumsum_num_turbines = cumsum(num_turbines))

# Column plot
ggplot(data = turbines_year,
       aes(x = year, y = cumsum_num_turbines)) +
  geom_col(fill = "darkblue") +
  labs(title = "Cumulative number of wind turbines installed in Rogaland",
       y = "",
       x = "") +
  theme_bw()

How would you improve these figures? For example, we can change the background colour, add subtitles and captions, change font sizes, etc. Can you generate your own figure? Why do you think it looks better? As an example I have created this one, what do you think?

caption_text <- "Data source: The Norwegian Water Resources and Energy Directorate (NVE)\nAuthor: Javier Elío (@Elio_javi) - Western Norway University of Applied Sciences"

# Column plot
ggplot(data = power_year,
       aes(x = year, y = cumsum_power_MW)) +
  geom_col(fill = "#0072B2") +
  labs(title = "Wind energy in Rogaland (Norway)",
       subtitle = "Cumulative power installed capacity in MW",
       caption = caption_text,
       y = "",
       x = "") +
  expand_limits(y = c(0, 1700)) +
  theme_bw() +
  theme(
    # Title and captions
    plot.title = element_text(size = 15, colour = "darkblue", face = "bold"),
    plot.subtitle = element_text(size = 10, colour = "grey25"),
    plot.caption = element_text(size = 10, colour = "grey25"),
    # Background colour
    plot.background = element_rect(fill = "linen", colour = NA),
    panel.background = element_rect(fill = "grey85", colour = NA)
  ) +
  # Add arrow
  annotate(
    'curve',
    x = as.Date("2016-01-01"), # Play around with the coordinates until you're satisfied
    y = 800,
    yend = 1600,
    xend = as.Date("2021-01-01"),
    linewidth = 1.5,
    curvature = 0.3,
    col = "#D55E00",
    arrow = arrow(length = unit(0.5, 'cm'))
  ) +
  # Add text
  annotate(
    'text',
    x = as.Date("2012-06-01"),
    y = 1300,
    label = "The installed capacity has nearly\ndoubled between 2017 and 2021.\nWhat will the limit be?",
    size = 3.5,
    hjust = 0
  )

Exercise 3

Objetives

• Get wind data from NORA3 link.

The R-script is based on matlab fucntions link

Solution

library(terra)
library(tidyterra)
library(eurostat)
library(sf)
library(tidyverse)
library(purrr)
library(patchwork)
library(ncdf4)

# # EU map
# box <- st_bbox(c(xmin = -20, xmax = 20, ymax = 45, ymin = 80),
#                 crs = st_crs(4326)) |> 
#   st_as_sfc() |> 
#   st_transform(3035)

eu_countries <- get_eurostat_geospatial(resolution = 10, 
                                        nuts_level = 0, 
                                        year = 2016,
                                        crs = "3035") 

Read data directly from the web (wiothout dowloading the data to a local folder

#' Get the data of wind speed and direction at specific heihts from a raster
#' @param .r NORA3 data (SpatRaster)
#' @param .height Height to get the data (20, 50, 100, 250, 750)

get_wind_height <- function(.r, .height = 100){
  
  # Velocity vector (y, y)
  ux = subset(.r, paste0("x_wind_", .height, "m"))
  uy = subset(.r, paste0("y_wind_", .height, "m"))
  
  # Calculate magnityd (mag) and direction (dir)
  u_mag = sqrt(ux^2 + uy^2)
  names(u_mag) <- "magnitude"
  u_dir = terra::atan2(y = uy, x = ux) * 180/pi
  names(u_dir) <- "direction"

  # Generate raster 
  u = c(u_mag, u_dir) 
  
  # Output as one raster
  return(u)
  
}

# Function for dowloading wind data from NORA3 
get_wind_z <- function(.year,
                       .month,
                       .day,
                       .hour_group,
                       .lead_time){
  
  # URL of the data
  nora3_url <- paste0("https://thredds.met.no/thredds/dodsC/nora3/",
                      .year,
                      "/",
                      .month,
                      "/",
                      .day, 
                      "/",
                      .hour_group,
                      "/fc", 
                      .year,
                      .month, 
                      .day,
                      .hour_group,
                      "_",
                      .lead_time,
                      "_fp.nc")
  
  # Open the netCDF file
  ncin <- ncdf4::nc_open(nora3_url)
  
  # Get coordinate  variables
  lon <- ncdf4::ncvar_get(ncin,"x")
  lat <- ncdf4::ncvar_get(ncin,"y")
  
  # Get time
  time <- ncdf4::ncvar_get(ncin,"time")
  
  # Get wind speed at 10 m above ground (height4)
  dname <- "wind_speed"
  ws10_array <- ncdf4::ncvar_get(ncin,dname)
  dlname <- ncdf4::ncatt_get(ncin,dname,"standard_name")
  dunits <- ncdf4::ncatt_get(ncin,dname,"units")
  fillvalue <- ncdf4::ncatt_get(ncin,dname,"_FillValue")
  # replace netCDF fill values with NA's
  ws10_array[ws10_array == fillvalue$value] <- NA
  
  # Get wind direction at 10 m above ground (height4)
  dname <- "wind_direction"
  wd10_array <- ncdf4::ncvar_get(ncin,dname)
  dlname <- ncdf4::ncatt_get(ncin,dname,"standard_name")
  dunits <- ncdf4::ncatt_get(ncin,dname,"units")
  fillvalue <- ncdf4::ncatt_get(ncin,dname,"_FillValue")
  # replace netCDF fill values with NA's
  wd10_array[wd10_array == fillvalue$value] <- NA
  
  # x_wind_z[x,y,height2,time] 
  # ncdf4::ncvar_get(ncin,"height2") -- 20  50 100 250 500 750 m above ground
  dname <- "x_wind_z"
  xh2_array <- ncdf4::ncvar_get(ncin,dname)
  dlname <- ncdf4::ncatt_get(ncin,dname,"standard_name")
  dunits <- ncdf4::ncatt_get(ncin,dname,"units")
  fillvalue <- ncdf4::ncatt_get(ncin,dname,"_FillValue")
  # replace netCDF fill values with NA's
  xh2_array[xh2_array == fillvalue$value] <- NA
  
  # y_wind_z[x,y,height2,time]
  dname <- "y_wind_z"
  yh2_array <- ncdf4::ncvar_get(ncin,dname)
  dlname <- ncdf4::ncatt_get(ncin,dname,"standard_name")
  dunits <- ncdf4::ncatt_get(ncin,dname,"units")
  fillvalue <- ncdf4::ncatt_get(ncin,dname,"_FillValue")
  # replace netCDF fill values with NA's
  yh2_array[yh2_array == fillvalue$value] <- NA
  
  
  # create dataframe with values
  df <- expand.grid(lon,lat) |> 
    as_tibble() |> 
    dplyr::rename_with(~ c("lon", "lat"), 1:2) |> 
    # Add wind speed and direction at 10 m
    dplyr::mutate(wind10_mag = as.vector(ws10_array),
                  wind10_dir = as.vector(wd10_array)) |> 
    # Add wind speed at h2 
    dplyr::mutate(
      # x_wind_z
      x_wind_20m  = as.vector(xh2_array[ , , 1]),
      x_wind_50m  = as.vector(xh2_array[ , , 2]),
      x_wind_100m = as.vector(xh2_array[ , , 3]),
      x_wind_250m = as.vector(xh2_array[ , , 4]),
      x_wind_500m = as.vector(xh2_array[ , , 5]),
      x_wind_750m = as.vector(xh2_array[ , , 6]),
      # y_wind_z 
      y_wind_20m  = as.vector(yh2_array[ , , 1]),
      y_wind_50m  = as.vector(yh2_array[ , , 2]),
      y_wind_100m = as.vector(yh2_array[ , , 3]),
      y_wind_250m = as.vector(yh2_array[ , , 4]),
      y_wind_500m = as.vector(yh2_array[ , , 5]),
      y_wind_750m = as.vector(yh2_array[ , , 6])
    )
  
  # Create Raster with all data
  r_crs <- "+proj=lcc +lat_0=66.3 +lon_0=-42 +lat_1=66.3 +lat_2=66.3 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs"
  r <- tidyterra::as_spatraster(df, crs = r_crs, digits = 4)
  # Add time 
  time(r) <- rep(lubridate::as_datetime(time, tz = "UTC"), times = length(names(r)))
  
  # get wind speed and direction at all heights
  ff <- function(.height) { get_wind_height(r, .height) }
  u_height_list <- c(20, 50, 100, 250, 750)
  u_height <- map(u_height_list, ff)
  names(u_height) <- c("wind20", "wind50", "wind100", "wind250", "wind750")
  u_height <- u_height |>
    rast() |> 
    rename_with( ~ gsub("_1", "_mag", .x, fixed = TRUE)) |> 
    rename_with( ~ gsub("_2", "_dir", .x, fixed = TRUE))
  
  # Retrurn SpatRaster object
  rr <- c(tidyterra::select(r, wind10_mag,  wind10_dir), u_height)
  
  return(rr)
  
} 

wind_nora3 <- get_wind_z(.year = "2018",
                         .month = "11",
                         .day = "08",
                         .hour_group = "00",
                         .lead_time = "004")

wind_nora3

class       : SpatRaster 
dimensions  : 1489, 889, 12  (nrow, ncol, nlyr)
resolution  : 2999.875, 2999.969  (x, y)
extent      : 776860.9, 3443750, -1271977, 3194976  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_0=66.3 +lon_0=-42 +lat_1=66.3 +lat_2=66.3 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs 
source(s)   : memory
names       :   wind10_mag,   wind10_dir,   wind20_mag, wind20_dir,   wind50_mag, wind50_dir, ... 
min values  :  0.004474435, 3.791680e-04,  0.003715257,  -179.9983,  0.006306755,  -179.9968, ... 
max values  : 26.079982758, 3.599998e+02, 27.383629015,   179.9978, 29.652804073,   179.9996, ... 
time        : 2018-11-08 04:00:00 UTC 

ggplot() +
  geom_spatraster(data =  select(wind_nora3 , ends_with("mag"))) +
  facet_wrap(~lyr, ncol = 2) +
  scale_fill_whitebox_c(name = "u [m/s]", palette = "viridi") +
  geom_sf(data = eu_countries, fill = NA, colour = "white") +
  coord_sf(expand = FALSE, crs = sf::st_crs(3035)) +
  scale_x_continuous(limits = c(3000000, 5000000)) +
  scale_y_continuous(limits = c(3600000, 5300000)) +
  labs(title = terra::time(wind_nora3)[1]) +
  theme_bw()

Figure 4

Extract values from one point (e.g., near Utsira - Coordinates from google maps [WGS84]: 59.346588 (lat), 4.899842 (long)).

# Point 
xy <- cbind(4.899842, 59.346588)
point <- vect(xy, crs="+proj=longlat +datum=WGS84") |> 
  project(crs(wind_nora3))

# Extract values
p_wind <- terra::extract(wind_nora3, point) |>  
  select(ends_with("_mag")) |> 
  pivot_longer(cols = everything()) |> 
  rename(height_m = name,
         wind_m_s = value) |> 
  mutate(height_m = readr::parse_number(height_m))
  

# Plot
ggplot() + 
  geom_point(data = p_wind,
             aes(x = wind_m_s, y = height_m ),
             col = "blue") +
  labs(x = "u [m/s]",
       y = "height [m]",
       title = "Vertical wind profile at point: 59.346588 (lat), 4.899842 (long)", 
       subtitle = terra::time(wind_nora3)[1]) +
  theme_bw()

Interpolation wind profile (Solbrekke and Sorteberg 2022):

• Exponential relation:

  \[ u_{z_2}(t) = u_{z_1}(t) \Big(\frac{z_2}{z_1}\Big)^{\alpha(t)} \]

• Exponential power law coefficient:

  \[ \alpha(t) = \frac{ln\Big( \frac{u_{z_2}(t)}{u_{z_1}(t)} \Big)}{ln\Big( \frac{z_2}{z_1} \Big)} \]

get_inter_wind_profile <- function(x){
  
  # Vector with alpha depending on height
  alpha <- rep(NA, 5)
  for(i in seq_along(alpha)) {
    alpha[i] = as.numeric( (log(x[i, 2]/x[i+1, 2]) / log(x[i, 1]/x[i+1, 1])) )
  }
  
  # Create data frame for interpolation (steps = 1) 
  wind_inter <- tibble(height_m = seq(10, 750, 1)) |> 
    mutate(wind_m_s = NA) |> 
    mutate(alpha = case_when(
      height_m <= 20 ~ alpha[1],
      height_m >  20 &  height_m <=  50 ~ alpha[2],
      height_m >  50 &  height_m <= 100 ~ alpha[3],
      height_m > 100 &  height_m <= 250 ~ alpha[4],
      height_m > 250 &  height_m <= 750 ~ alpha[5]
    ))  
  
  # Interpolation (based on u [m/s] at 10 m)
  wind_inter$wind_m_s[1] = x[x["height_m"] == 10, ]$wind_m_s
  for(i in 2:length(wind_inter$height_m)) {
    wind_inter$wind_m_s[i] = wind_inter$wind_m_s[i-1] * (wind_inter$height_m[i] / wind_inter$height_m[i-1])^wind_inter$alpha[i]
  }
  
  return(wind_inter)
  
}
  
wind_inter <- get_inter_wind_profile(p_wind)

# Plot
ggplot() + 
  geom_point(data = p_wind,
             aes(x = wind_m_s, y = height_m ),
             col = "blue") +
  geom_line(data = wind_inter,
             aes(x = wind_m_s, y = height_m ),
             col = "red") + 
  labs(x = "u [m/s]",
       y = "height [m]",
       title = "Vertical wind profile at point: 59.346588 (lat), 4.899842 (long)", 
       subtitle = terra::time(wind_nora3)[1]) +
  theme_bw()

Exercise 4

Objetives

Time series wind using NORA3 data

Analyse wind speed and direction over time in a point at 100 m hight.

Solution

library(terra)
library(tidyterra)
library(eurostat)
library(sf)
library(tidyverse)
library(lubridate)
library(purrr)
library(patchwork)
library(ncdf4)
library(gt)

Read data that have been previously downloaded into a local file.

wind_100m <- readRDS("data/wind_100m.rds")

Hourly analysis

Use average hourly data directly from NORA3.

ggplot(data = wind_100m, aes(x = time, y = u_mag)) +
  geom_point() +
  # geom_smooth(span = 0.3, se = TRUE) +
  # Cut-in
  geom_hline(yintercept = 3.5, col = "#009E73", linetype = "dashed") + 
  geom_text(x = as_datetime("2018-01-01 04:00:00"),
            y = 4,
            label = "Cut-in",
            colour = "#009E73") +
  # Cut-off
  geom_hline(yintercept = 25, col = "#D55E00", linetype = "dashed") +
  geom_text(x = as_datetime("2018-01-01 04:00:00"),
            y = 25.5,
            label = "Cut-off",
            colour = "#D55E00") +
  # Format labels and title
  labs(x = "Time",
       y = "u [m/s]",
       title = "Hourly wind speed at 100 m height", 
       subtitle = "Coordinates (WGS84): Long = 4.899842, lat = 59.346588",
       caption = "Data: NORA3") +
  coord_cartesian(ylim = c(0,30)) +
  theme_bw()

Figure 5: Hourly wind speed (m/s)

Histogram

# Plot histogram 
ggplot() +
  geom_histogram(data = wind_100m, aes(x = u_mag, y =..density..),
                 binwidth = 0.5, 
                 fill = "#D55E00", 
                 colour = "grey") +
  # Format labels and title
  labs(x = "u [m/s]",
       title = "Histogram of hourly wind speed at 100 m height", 
       subtitle = "Point coordinates (WGS84): Long = 4.899842, lat = 59.346588",
       caption = "Data: NORA3") +
  coord_cartesian(xlim = c(0, 65)) +
  theme_bw() 

Figure 6: Wind speed histogram

Fit a Weibull distribution

# Fit weibul distribution
fit <- fitdistrplus::fitdist(wind_100m$u_mag,"weibull")

# Put parameter in a table (inset element)
weibull <- tibble(param = rownames(as.data.frame(fit$estimate)),
          estimate = fit$estimate,
          std_Error = fit$sd
          )
# Generate values for addind the density curve
fit_values = rweibull(1000000, shape = fit$estimate[1], scale = fit$estimate[2]) |> 
  as_tibble()

# Plot histogram together with the weibull distribution 
ggplot() +
  geom_histogram(data = wind_100m, aes(x = u_mag, y =..density..),
                 binwidth = 0.5, 
                 fill = "#D55E00", 
                 colour = "grey") +
  # Addd Weibull
  geom_density(data = fit_values,
               aes(x = value, y = ..density..),
               colour = "#0072B2",
               linewidth = 1) +
  # Format labels and title
  labs(x = "u [m/s]",
       title = "Weibull distribution of hourly wind speed at 100 m height", 
       subtitle = "Point coordinates (WGS84): Long = 4.899842, lat = 59.346588",
       caption = "Data: NORA3") +
  coord_cartesian(xlim = c(0, 65)) +
  theme_bw() +
  # Add inset table
  inset_element(gt(weibull), 0.60, 0.75, 0.95, 0.95)

Figure 7: Weibull distribution

Time bellow/above thresholds:

# Time in production 
wind_100m <- wind_100m |> 
  select(time, u_mag, u_dir) |> 
  mutate(production = case_when(
    u_mag < 3.5 ~ FALSE,
    u_mag >  25 ~ FALSE,
    .default = TRUE)
    ) 

wind_100m |> 
  summarise(
    total_hours = n(),
    total_prod = sum(production),
    perc = 100 * mean(production)
  ) |> 
  gt() |> 
  fmt_number(columns = perc,
             decimals = 1)

Table 1: Percentage of time that a wind turbine is producing energy based on average hourly wind speed data (cut-in = 3.5 m/s, and cut-off = 25 m/s)

total_hours	total_prod	perc
744	688	92.5

Ten minutes resolution

Create simulated data (n = 1000) based on NORA3 hourly data. Simulated data follow a Weibull distribution with the same global shape but mean values of the hourly data.

set.seed(36)
n_sim = 1000


# Shape
wd_shape <- fit$estimate[1]

# Initial data frame  
sim_wind_100m <- tibble(time = seq(ymd_hm("2018-01-01 04:00"),
                                   ymd_hm("2018-02-01 03:50"), 
                                   by = "10 mins")) |> 
  # Add mean values and simulate weibull 
  left_join(wind_100m) |> 
  dplyr::select(time, u_mag) |> 
  fill(u_mag)


# Simulations
wd_shape <- fit$estimate[1]
sim <- matrix(NA, length(sim_wind_100m$u_mag), n_sim)
u_mag = sim_wind_100m$u_mag
for(i in 1:n_sim) {
  for(j in seq_along(u_mag)) {
    sim[j, i] = rweibull(1, shape = wd_shape, scale = u_mag[j] / gamma(1+(1/wd_shape)))
  }
}

sim <- sim |> 
  as_tibble() |> 
  rename_with(~ paste0("sim_", .)) |> 
  rename_with(~ gsub("V", "", .x, fixed = TRUE))
  
# Add to data frame
sim_wind_100m <- sim_wind_100m |> 
  bind_cols(sim)

# Plot example of 4 simulations
sim_wind_100m |> 
  select("sim_1", "sim_2", "sim_3", "sim_4") |> 
  pivot_longer(cols = starts_with("sim"))  |> 
  ggplot(aes(x = value)) +
  geom_histogram(aes(y =..density..),
                 bins = 50,
                 fill = "#D55E00",
                 colour = "grey") +
  geom_density() +
  labs(title = "Wind distribution in simulated data",
       subtitle = "Example of the first four simulations",
       x = "u [m/s]") +
  facet_wrap(~ name, ncol = 2) +
  theme_bw() 

Figure 8

# Exampla first four simulations
sim_wind_100m |> 
  select("time","u_mag", "sim_1", "sim_2", "sim_3", "sim_4") |> 
  pivot_longer(cols = starts_with("sim")) |> 
  ggplot() +
  geom_point(aes(x = time, y = value),
             colour = "#0072B2",
             size = 0.6,
             alpha = 0.5) + 
  # Hourly points
  geom_line(data = wind_100m, aes(x = time, y = u_mag), col = "red") +
  # Cut-in
  geom_hline(yintercept = 3.5, col = "#009E73", linetype = "dashed") +
  # Cut-off
  geom_hline(yintercept = 25, col = "#D55E00", linetype = "dashed") +
  # Format labels and title
  labs(x = "Time",
       y = "u [m/s]",
       title = "Simulated wind speed at 100 m height (10 min resolution)", 
       subtitle = "Exampla of the first four simulations",
       caption = "Red line indicates hourly mean value from NORA3") +
  coord_cartesian(ylim = c(0, 90)) +
  facet_wrap(~ name, ncol = 2) +
  theme_bw()

Figure 9: Hourly wind speed (m/s) of simntetic data

sim_wind_100m |> 
  pivot_longer(cols = starts_with("sim")) |> 
  ggplot() +
  # geom_point(aes(x = time, y = value),
  #            colour = "#0072B2",
  #            size = 0.01,
  #            alpha = 0.05) + 
  geom_hex(aes(x = time, y = value), bins = 500) +
  viridis:: scale_fill_viridis(trans = "log10") +
  # Hourly points
  geom_line(data = wind_100m,
            aes(x = time, y = u_mag), 
            col = "red",
            alpha = 0.5) +
  # Cut-in
  geom_hline(yintercept = 3.5, col = "#009E73", linetype = "dashed") +
  # Cut-off
  geom_hline(yintercept = 25, col = "#D55E00", linetype = "dashed") +
  # Format labels and title
  labs(x = "Time",
       y = "u [m/s]",
       title = "Simulated wind speed at 100 m height (10 min resolution)", 
       subtitle = "All simulations",
       caption = "Red line indicates hourly mean value from NORA3") +
  coord_cartesian(ylim = c(0, 90)) +
  theme_bw()

Hourly average vs. NORA3

sim_wind_100m |>  
  mutate(date_hour = format(time, "%Y-%m-%d %H")) |> 
  relocate(date_hour, .after = time) |> 
  pivot_longer(cols = starts_with("sim")) |> 
  group_by(date_hour) |> 
  summarise(
    n = n(),
    time = first(time),
    u_mag = first(u_mag),
    sim_mean = mean(value)
    ) |> 
  ungroup() |> 
  ggplot() +
  geom_point(aes(x = u_mag, y = sim_mean), alpha = 0.5) +
  labs(title = "Hourly wind speed [m/s]",
       y = "Simulated data",
       x = "NORA3 data") +
  geom_abline(intercept = 0, slope = 1, col = "red") +
  theme_bw()

Figure 10: Relationship between the hourly wind speed calculated in the simulated data and NORA3 data values

Percentage of time that a wind generator is producing energy in the 10 minute simulated data:

prod_10min <- sim_wind_100m |> 
  pivot_longer(cols = starts_with("sim")) |> 
  group_by(name) |>
  mutate(production = case_when(
    value < 3.5 ~ FALSE,
    value >  25 ~ FALSE,
    .default = TRUE)
    ) |> 
  summarise(
    total_hours = n() * 10/60,
    total_prod = sum(production) * 10/60,
    perc = 100 * mean(production)
    ) |> 
  ungroup()


# Plot
prod_10min |> 
  ggplot(aes(x = perc)) +
  geom_histogram(aes(y =..density..),
                 bins = 50,
                 fill = "#D55E00",
                 colour = "grey") +
  geom_density() +
  labs(x = "Perc. [%]",
       title = "Hystogram of the percentage of hours in production",
       subtitle = "Wind speed between cut-in = 3.5 m/s and cut-off = 25 m/s") +
  theme_bw() 

Figure 11: Percentage of time that a wind generator is producing energy based on simulated data (cut-in = 3 m/s and cut-off = 25 m/s)

prod_10min |> 
  select("total_hours", "total_prod", "perc") |> 
  purrr::map_dfr(.f = summary, .id = "value") |> 
  mutate_at(vars(-value), as.numeric) |> 
  gt() |> 
  fmt_number(columns = vars(-value),
             rows = c(1,2),
             decimals = 0) |> 
  fmt_number(columns = vars(-value),
             rows = 3,
             decimals = 1)

Table 2: Percentage of time that a wind generator is producing energy based on simulated data (cut-in = 3 m/s and cut-off = 25 m/s)

value	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
total_hours	744	744	744	744	744	744
total_prod	554	565	568	568	570	581
perc	74.5	75.9	76.3	76.3	76.7	78.1

Analyse variability in the data.

Solbrekke, Ida Marie, and Asgeir Sorteberg. 2022. “NORA3-WP: A High-Resolution Offshore Wind Power Dataset for the Baltic, North, Norwegian, and Barents Seas.” Scientific Data 9 (1): 362. https://doi.org/10.1038/s41597-022-01451-x.